Dynamical bifurcation as a semiclassical counterpart of a quantum phase transition 
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We illustrate how dynamical transitions in nonlinear semiclassical models can be recognized as phase 
transitions in the corresponding - inherently linear - quantum model, where, in a Statistical Mechanics 
framework, the thermodynamic limit is realized by letting the particle population go to infinity at fixed 
size. We focus on lattice bosons described by the Bose-Hubbard (BH) model and Discrete Self-Trapping 
(DST) equations at the quantum and semiclassical level, respectively. After showing that the gaussianity 
of the quantum ground states is broken at the phase transition, we evaluate finite populations effects 
introducing a suitable scaling hypothesis; we work out the exact value of the critical exponents and we 
provide numerical evidences confirming our hypothesis. Our analytical results rely on a general scheme 
obtained from a large-population expansion of the eigenvalue equation of the BH model. In this approach 
the DST equations resurface as solutions of the zeroth-order problem. 



The emergence of nonlinear phenomena typical of semi- 
classical models from more fundamental, inherently linear 
quantum models represents a basic question which has at- 
tracted a fair amount of attention in the literature. The large 
number of experiments with ultracold atoms trapped in op- 
tical lattices evidencing both quantum OJ-Q] and nonlin- 
ear |H-i|] effects make the the Bose-Hubbard model lE^-fllll 
an important playground in this respect. A general form 
for the Bose-Hubbard Hamiltonian, describing interacting 
bosons hopping on the sites of a discrete structure, is 
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where j is a site label, Ajt is the adjacency matrix of 
the structure, dj is an on-site (annihilation) boson opera- 
tor, [cij, a\\ = 5j£. As to the parameters, U > is the 
strength of the on-site interaction, whose attractive or re- 
pulsive character is dictated by s = ±1, and J is the hop- 
ping amplitude. The semiclassical counterpart of Eq. £[]) 
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can be obtained by trading the quantum operator dj for a C- 
number a>j, whose square modulus and angle represent the 
boson population and macroscopic phase attached to site 
j, respectively. Actually, the complex variables oij gov- 
erned by Hamiltonian d2]l are the variational para meters of 
a coherent-state ansatz for Hamiltonian (OQ) fl2T4l4ll . The 
nonlinear equation of motion ensuing from Hamiltonian 
3 are known as discrete self-trapping (DST) equations 
1511 . and the relevant normal modes are often compared 
to the eigenstates of Hamiltonian (TJ. 

One of the most striking features stemming from the 
nonlinear character of DST equations is the occurrence 
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tty II 1511 : for (attractive) interactions among bosons ex- 
ceeding a critical value, the uniform solution of the DST 
equations on translation invariant lattices becomes unsta- 
ble. This critical value is expected to coincide with the 
threshold for spatial localization in the ground-state of the 
system (i.e. soliton formation). The quantum counterpart 
of this well known semiclassical feature has been investi- 
gated — mostly on two- and three-site lattices at finite pop- 
ulation — by employing several indicators, such as energy 
gaps lfl6h . number fluctuations [17-20]] condensate frac- 
tion IU4II . occupation probability in the Fock space JI21I1. lo - 
calization width ll22ll . fidelit y C3ll . entanglement I24l427ll . 
Bethe ansatz techniques Il24ll28ll . Here we show that this 
semiclassical dynamical transition can be seen as a gen- 
uine quantum phase transition - characterized by a van- 
ishing gap - in which the thermodynamic limit is realized 
at fixed lattice size by letting the bosonic population go 
to infinity. We connect the order of the phase transitions 
to the bifurcation pattern characterizing the solutions to 
the DST equations. Specifically, we find that on ID lat- 
tices comprising L = 2 and L > 6 sites the modulational 
instability threshold coincides with a second-order critical 
point, while for L = 3, 4, 5 only first-order transitions are 
present. In the second-order case, we determine the critical 
exponents characterizing the divergence of the fluctuations. 
Also, we analyze the crossover heralding the phase transi- 
tion at finite population by verifying the finite-size scaling 
hypothesis typical of Statistical Mechanics. In the same 
framework we show that at finite population the spectrum 
is not gapless at the transition point, but exhibits the ex- 
pected avoided level crossing pattern. 

In carrying out our analysis we develop a fully analyti- 
cal general scheme based on a large-population expansion 
of the secular equation for the BH Hamiltonian on a L- 
site lattice. This expansion results in a Schrodinger-like 
equation (SLe) set in a L-dimensional space directly re- 
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lated to the Fock space of the original quantum problem. In 
this framework the DST equations emerge naturally as the 
equations for the extrema of the potential part of the SLe, at 
the lowest order of the expansion. The analysis of the small 
oscillations about the local minima of such potential gives 
access to an entire class of solutions of the original quan- 
tum problem, which improve significantly on the coherent 
states employed as trial wave functions in the variational 
derivation of the DST equations. Similar to the coherent 
state, away from criticality our solution is a Gaussian func- 
tion centered at a particular occupation-number Fock state 
dictated by the time-independent DST equations. However, 
our perturbative approach provides an exact calculation of 
the Gaussians widths. This gives us access to the correct 
quantum fluctuations in the system. Furthermore, our pic- 
ture captures situations where the solutions of the SLe are 
not Gaussian, which correspond to the above mentioned 
quantum critical points. 

Our analysis starts from the eigenvalue equation for the 
BH model, H\W) = E\^f), where the eigenstate |*) has 
been expanded over the Fock space of occupation numbers, 

|*>=X>(aO|aO, I£>=n-7^rl°>- (3) 

x 3=1 K lyX j)- 

The labels of the Fock state \x) have been conveniently 
normalized, so that Xj € [0, 1] independent of the to- 
tal boson population N, which is conserved owing to the 
commutation relation [H, J2j a ] a j] = 0- The prime on 
the summation symbol in the first of Eqs. (0) signals that 
Xj = kjN -1 , with kj £ No, and Y^j x j = 1- When recast 
in terms of the expansion coefficients ip(x), the eigenvalue 
equation can be seen as a discrete equation on the vertices 
of a mesh grid covering the portion of (L — l)-dimensional 
hyperplane defined by J2j x j = 1 an ^ %j & [0, 1]. The 
N — > oo limit plays the twofold role of a thermody- 
namic limit for the size of the fixed-number Fock space, 
and of a continuous limit for the mesh grid. Assum- 
ing that the expansion coefficients can be seen as a con- 
tinuous function, V( x i; x 2, • • • , Xi + -/V -1 , • • • ,x N ) — 
tp(xi, X2, ■ ■ ■ ,X£,-- - ,x N ) = 0(N^ 1 ) the eigenvalue 
equation becomes a SLe of the form \U — T>]tp(x) = 
Eip(x) where, to the leading order in iV -1 , the self-adjoint 
operators are 

Uix) = ^2 ( sx ) + VjXj) - 2t ^2 Aji^/xjXi (4) 

j it 

t x - d d 

T^\ x ) = ~zt^ 2_^, i iT )i \V x l^TT f r (5) 

it ~ Xl > ^ X 'J ~ Xt > 

The effective parameters appearing in Eqs. (01) and (O 
are r = J/UN, v s = v 3 /UN and E = E/UN 2 . 
We remark that a similar approach has been adopted in 
Refs. Hzl [01 IH |H, HI, mostly for the two-site case. 
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FIG. 1: (color online) A. density plot for \ip(xi, X2, l—x% — X2)\ 2 
for a mid-spectrum eigenstate of a L — 3 lattice containing 
iV = 600 bosons, at r = 0.26 > t\ oc . The discontinuous char- 
acter of the expansion coefficients is apparent; B. Normalized 
number fluctuations vs. effective hopping amplitude on a lattice 
comprising L = 6 sites for increasing N. The thick dashed and 
dotted lines are the theoretical prediction obtained from Eq. ^ 
and from the coherent state approach fl2T - ll4ll . 

Note that the leading terms in D(x) and U(x) are of or- 
der iV° and iV -2 , respectively. This suggests that signif- 
icant solutions to the above differential equation might be 
strongly localized in the vicinity of the local minima of 
tl{x). Thus, a satisfactory description of a subset of eiges- 
tates of Hamiltonian £Q can be obtained through the anal- 
ysis of the "small oscillations" about such minima. 

The stationarization of the leading term of Eq. I© results 
in the set of equations 

sxj y/x~ + vj y/x~ - r Aji y/x~t = X^/x~, (6) 

i 

where A is a Lagrange multiplier enforcing the constraint 
on x. Note that, upon setting ^/xj = ctj, Eqs. (O are 
equivalent to the time-independent DST equations, which 
thus emerge in our description without any variational 
ansatz. 

If the above stationary points correspond to non- 
degenerate locally quadratic minima of U, the problem is 
mapped onto a quantum harmonic oscillator. In particu- 
lar, the expansion coefficients ip(x) of the (local) ground 
state is a Gaussian whose width is in principle accessible. 
As we mention the trial-wavefunction in the coherent-state 
approach of Refs. lfl3ll22ll is also a Gaussian, whose center 
is the Fock state corresponding to the relevant stationary 
solution of the DST equations. However, once the cen- 
ter is chosen, the width of such Gaussian is constrained 
by the coherent- state structure. Conversely, in our picture, 
the analysis of the local harmonic potential provides the 
correct width for the Gaussian function, which results the 
exact description of the scaling of quantum fluctuations in 
the N — > oo limit. Also, while coherent trial wavefunc- 
tions are inherently Gaussian, our description allows for 
non-Gaussian eigenstates. It should be remarked that while 
our approach is surely applicable to low-lying eigenfunc- 
tion of Hamiltonian (QQ|, the assumption of continuity may 
not hold for higher energies, as exemplified in Fig.[TJ 
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We now specialize the above discussion to the study of 
the localization/delocalization transition for attractive bo- 
son on ID, translation-invariant (Vj = 0) lattice com- 
prising L sites. As it is well known, the uniform state 
Xj = L _1 is always a solution to Eqs j6), and it coincides 
with the absolute minimum ofU(x) l22ll3~lll for attractive 
interactions and r > r loc as well as for repulsive interac- 
tions at any r. Below such localization threshold the min- 
imum of U(x) is L— fold degenerate and the correspond- 
ing solutions of Eqs. (O spontaneously break the transla- 
tional symmetry due to the nonlinear interaction term. For 
L > 5 the localization threshold coincides with the criti- 
cal value below which the uniform state becomes modula- 
tionally unstable, Ti oc = r m .i. = [2L sin 2 (7r/L)]~ 1 . The 
same happens in the two-site case. For L = 3, 4, 5 the low 
lying stationary points ofV((x) exhibit a more complex bi- 
furcation pattern, and there is a region of metastability of 
the uniform state, r m i . < T\ oc . As a consequence the L 
equivalent, symmetry breaking local minima of U(x) do 

symmetri 

of the ground state of Eq. (fl]l is particularly simple in the 
region r > t 1oc , where the absolute minimum of lA(x) 
is always locally harmonic. In order to analyze the small 
oscillations about such minimum we introduce the small 
deviations ij k 



not merge continuously with the symmetric minimum at 
t = r ioc> but disappear abruptly 11221 l3~ill . The analysis 
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This choice decouples the leading terms of the equation 
\U — T>\ip = Etp into L—l independent equations, 
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where if, = U k fa, E = 4r/(LA) 2 £ fc S k /r k , and 
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In particular, the lowest energy gap of the original problem, 
Eq. £Q, is AE = AJ/L^J(sti + t)/(ttJV We note that 
the above results were obtained in Ref. ltl8ll within a Bo- 
goliubov approach and in Ref. ll29ll . where Eq. dU) was de- 
rived introducing the small deviations from the uniform so- 
lution, Eq. (Q, directly into the structure of the Fock states. 

Our method is however more general, since it can be 
in principle applied to any stable solution of the time- 
independent DST equation, emerging as a non-degenerate 
locally harmonic minimum of U(x). Again, the eigen- 
state of Hamiltonian dTJ is a Gaussian state whose square 
width vanishes as A -1 . Note that, unlike Eq. ©, the 
square width of the uniform solution in the coherent-state 
approach |l2- 14| is (NL)^ 1 , irrespective of r and k. 

The above calculation allows us to evaluate, the (normal- 
ized) number fluctuations A(r, N) = N((x 2 ) — (xj) 2 ) = 
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FIG. 2: (color online) Ground state for N = 1200 attractive 
bosons on a lattice comprising L = 3 sites, at r = 0.21 < t\ oc — 
0.25. A. Density plot of \%l}{x\,X2, 1 — x\ — X2)\ 2 showing that 
that it consists of the superposition of three peaks. B. Plot along 
the cut x\ = X2 (dashed white line in panel A.) showing that the 
peaks (black thin line) are well approximated by a Gaussian func- 
tion (green thick line). The discrepancy at the tails is a finite-size 
(population) effect. See Ref. l22ll for more plots. 



n E kf0 (vl) = n E M0 °K T \ where the factor N en " 

sures that this quantity remains finite in the large— A limit. 
Panel B in Fig. Q] shows that numerical estimates of A ob- 
tained from population quantum Monte Carlo simulations 
l32ll exhibit a very satisfactory agreement with our analyt- 
ical prediction in the large- A limit. 

In the study of the ground state for attractive interactions 
(s = —1), two situations do not fit in the above general 
picture, and require a separate discussion. For r < T\ oc 
the potential has L degenerate (harmonic) absolute minima 
while for r = r loc = r m ; , on lattices comprising L = 2 
or L > 5 sites, the (non degenerate) absolute minimum of 
the potential is not harmonic. In the first case, the expan- 
sion coefficients ip(x) are the superposition of L Gaussian 
functions having the same width and centered at the abso- 
lute minimum points. Fig. |2] illustrates an example of this 
situation for L = 3. Thus the r < r loc interval is anal- 
ogous to a first-order phase transition line, where several 
different phases coexist. 

Likewise the r = r loc = r m .i. case is analogous to a 
second-order critical point. Note indeed that, since r±i = 
r m i , the fluctuations of the first "small oscillation" modes 



diverge as ct ±1 (t») ~ r, 
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where 

1/2 



r — r m .i.. Also 



the gap vanishes as AE ~ t* 7 ^. This signals that the har- 
monic term of IA vanishes for r = r m .i l33ll , so that the 
remaining leading term is quartic, and one has [r^T^ — 
N~ 2 d 2 1 ' drj 2 _ l \ r ij)±i = E±iip±i. Introducing the rescaled 
variables fj = N 1/3 rj ±1 and £ = N 4 / 3 E ±1 we find that 
V'±i( r ?±i) = 9 7 (-^ 1 ^ 3? 7±i) ; where <p is the (non-Gaussian) 
ground state solution of [r ±1 7) 4 — d 2 /dfj 2 ]ip = £(p. This 
means that, for r* = at finite A, the system is not gapless 
as suggested by the Gaussian approximation, but exhibits 
the expected avoided level crossing. Straightforward cal- 
culations allow us to determine the critical exponent gov- 
erning the vanishing of the energy gap with increasing pop- 
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FIG. 3: (color online) A. Scaling of N(rjj.) for L — 6 and sev- 
eral TV. The numeric data points (symbols) are compared to the 
theoretical predictions (solid lines) in Eqs. ^ and d 1 0b . B. Data 
collapse of finite— N numerical data showing the correctness of 
the scaling hypothesis for the fluctuation of the first modes. 



ulation, AE oc JN 1 ' 3 . Likewise, we obtain 
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N' 1 for Ifcl > 1, (10) 



i.e. the singular part of the number fluctuation A is dic- 
tated by (77^), where the scaling of the \ k\ > 1 modes is 
the same as far from criticality. Panel A. in Fig.[3]illustrates 
the correctness of Eqs. ( PTOl) in the L = 6 case. Assuming 
that (77^) obeys the standard finite-size scaling relation, 

(rf ±1 ) = N~> /v f(N 1/v T*), from Eqs. © and QOj one eas- 
ily works out the critical exponents 7 = — 1 and v = 3/2. 
As we mention, here the boson population plays the role 
reserved for the system size in the usual scaling approach. 
Panel B. of Fig. [3] illustrates how numerical data for a lat- 
tice comprising L = 6 sites nicely collapse according to 
the above scaling hypothesis. We conclude by remarking 
that for 2 < L < 6, where r m j < r loc , there is no 2nd or- 
der critical point. The situation is similar to ferromagnetic 
Ising models with p-body interaction (p > 2) ll34Tl . 

In summary, we develop a general technique for the 
study of bosonic quantum system in the large occupation 
limit. In this framework the corresponding semiclassical 
model is recovered without the need of an ansatz for the 
structure of the quantum state, and the dynamical transi- 
tions thereof incontrovertibly emerge as genuine quantum 
phase transitions. We study in detail the localization tran- 
sition of the paradigmatic Bose-Hubbard model with at- 
tractive interactions, for which we calculate explicitly some 
significant critical exponents and perform a finite-size scal- 
ing analysis. Our findings are supported by extensive nu- 
merical simulations. 
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